StatQuest生物统计学专题 - 随机森林(2) R实例

下载数据数据清洗构建随机森林随机森林评价

本文将在R中进行随机森林的实际操作。

原始数据是使用年龄、性别、胸痛程度、血压等共14个参数来预测一个人是否罹患心脏病。原始数据位于http://archive.ics.uci.edu/ml/datasets/Heart+Disease。

用到的R包为ggplot2、randomForest。

library(ggplot2)
library(randomForest)

下载数据

使用read.csv函数读取网络数据。

url <- "http://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.cleveland.data"

data <- read.csv(url, header=FALSE)

数据清洗

修改列名

# 取回的数据没有列名
head(data)
# 手动修改列名
colnames(data) <- c(
  "age",  "sex",  "cp",   "trestbps", 
  "chol",   "fbs",    "restecg", 
  "thalach",   "exang",    "oldpeak", 
  "slope",  "ca",   "thal",  "hd" 
  ) # 共14个

head（data）

修改数据格式

从官网可以知道14个数据变量的含义，他们分别是：

  age
  sex: 0 = female, 1 = male
  cp: chest pain
          # 1 = typical angina,
          # 2 = atypical angina,
          # 3 = non-anginal pain,
          # 4 = asymptomatic
  trestbps: resting blood pressure (in mm Hg)
  chol: serum cholestoral in mg/dl
  fbs: fasting blood sugar greater than 120 mg/dl, 1 = TRUE, 0 = FALSE
  restecg: resting electrocardiographic results
          # 1 = normal
          # 2 = having ST-T wave abnormality
          # 3 = showing probable or definite left ventricular hypertrophy
  thalach: maximum heart rate achieved
  exang: exercise induced angina, 1 = yes, 0 = no
  oldpeak: ST depression induced by exercise relative to rest
  slope: the slope of the peak exercise ST segment
          # 1 = upsloping
          # 2 = flat
          # 3 = downsloping
  ca: number of major vessels (0-3) colored by fluoroscopy
  thal： this is short of thalium heart scan
          # 3 = normal (no cold spots)
          # 6 = fixed defect (cold spots during rest and exercise)
          # 7 = reversible defect (when cold spots only appear during exercise)
  hd： (the predicted attribute) - diagnosis of heart disease
          # 0 if less than or equal to 50% diameter narrowing
          # 1 if greater than 50% diameter narrowing

然后我们使用str(data)看一下此时的数据结构可以发现不规范的变量有，

• sex：应该是因子变量，0代表女，1代表男；

• cp、fbs、restecg、exang、slope：应该是因子变量；

• ca：应该使用“NA”表示缺失值，而不是“？”；

• thal：和ca情况一样；

• hd：因子变量，0是健康，1是心脏病；

按照上述情况，进行以下修改：

#先将“？”更换为“NA”
data[data=="?"]<-NA

#将sex修改为正确的因子变量，并处理好性别
data$sex <- factor(data$sex,levels = c(0,1), labels = c("F","M"))
#将cp、fbs、restecg、exang、slope修改为因子变量
data$cp <- as.factor(data$cp)
data$fbs <- as.factor(data$fbs)
data$restecg <- as.factor(data$restecg)
data$exang <- as.factor(data$exang)
data$slope <- as.factor(data$slope)

#由于ca刚开始使用的’？‘代表的缺失值，因此ca里面的变量都是字符类型的因子变量
#因此需要先转换为数值型，再进行因子变量转换
data$ca <- as.integer(data$ca)
data$ca <- as.factor(data$ca)
#同理，将thal做同样处理
data$thal <- as.integer(data$thal) 
data$thal <- as.factor(data$thal)

#将hd修改为因子变量，0为健康，1为不健康
data$hd <- ifelse(data$hd == 0, "Healthy", "Unhealthy")
data$hd <- as.factor(data$hd) 

构建随机森林

填补缺失值

rfImpute函数用于填补缺失值，随机森林的缺失值填补是根据相似度进行填补的一种迭代算法。

简单来说，缺失值会首先根据数据类别不同进行填充，也就是数值数据填充当前变量的中位数，分类数据填充当前变量的众数。然后构建随机森林，并根据随机森林来决定所有的记录之间的相似度，最后根据相似度和当前变量的其他数据对缺失值进行填补。

一般上述过程会迭代4-6次，4-6次被认为是比较合理的迭代次数。

R中ifImpute函数默认创建300棵决策树的随机森林。

最佳子集数目根据数据类别不同进行设定，数值数据为总变量数除以3，分类数据为总变量数的平方根。

# rfImpute填补缺失值
# hd~. 其中hd为相应变量，.为其他所有变量，其意义为使用其他所有变量预测hd
data.imputed <- rfImpute(hd ~ ., data = data, iter=6)

结果会输出每次迭代后的OOB值，越低越好。

构建随机森林

# 构建随机森林
# 默认创建500棵决策树的随机森林。
# 最佳子集数目根据数据类别不同进行设定，数值数据为总变量数除以3，分类数据为总变量数的平方根。

model <- randomForest(hd ~ ., data=data.imputed, proximity=TRUE)
# proximity参数不是必须的，加上后，则会输出proximity矩阵，此矩阵可用于热图或MDS（PCoA）

随机森林评价

决策树的数量

默认是创建500棵决策树，此时的OOB（out of bag）值可以用于评价随机森林的模型如何。

如果你不知道什么是OOB，请参阅上一节内容随机森林(1)。

我们可以看看此时从第1棵树到第500棵决策树时，OOB的变化趋势：

# model下err.rate中是OOB的数据，它有三列，分别是总OOB、健康人的OOB以及不健康人的OOB
# 创建数据框用于ggplot绘图
oob.error.data <- data.frame(
  Trees=rep(1:nrow(model$err.rate), times=3),
  Type=rep(c("OOB", "Healthy", "Unhealthy"), each=nrow(model$err.rate)),
  Error=c(model$err.rate[,"OOB"],
    model$err.rate[,"Healthy"],
    model$err.rate[,"Unhealthy"]))
# 绘图
ggplot(oob.error.data, aes(x=Trees, y=Error)) + geom_line(aes(color=Type))

可以看出，大概从150以后的OOB的值趋于稳定了，默认的500是非常稳健的数值了。

StatQuest-XXIV-1

最佳子集数目

默认子集数目为总变量数的平方跟，也就是13的平方根，约为3.6，所以默认的子集数目为3。

我们可以改变不同的子集数目以确认最佳子集数目是多少，比如可以看一下子集数目分别为1-10时的结果：

oob.values <- vector(length=10)
for (i in 1:10){
    temp.model <- randomForest(hd~.,data = data.imputed,mtry=i)
    oob.values[i] <- temp.model$err.rate[nrow(temp.model$err.rate),1]
}
oob.values
# Result:0.1749175 0.1716172 0.1716172 0.1815182 0.1749175 0.1848185 0.1881188 0.1947195 0.1914191 0.1848185

可以发现最低的OOB确实是子集数目为3。

专题以往文章

参考资料

StatQuest课程：https://statquest.org/video-index/

---

猜你喜欢

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。

▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。